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ABSTRACT 

An approach to calculating optimal, gliding 
flight paths of the type associated with the space 
shuttle's transition from entry to cruising flight is 
presented* Kinetic energy and total energy (per 
unitweight) replace velocity and time in the dynamic 
equations, reducing the dimension and complexity 
of the problem* The capability for treating integral 
and terminal penalties (as well as Mach number 
effects) is retained in the numerical optimization; 
hence, stability and control boundaries can be ob- 
served as trajectories to the desired final energy, 
flight path angle, and range are determined* Numer- 
ical results show that the "jump' 1 to the "front-side 
of the L ID curve" need not be made until the end 
of the transition and that the dynamic model pro- 
vides a conservative range estimate* Alternatives 
for real-time trajectory control are discussed. 


INTRODUCTION 

The space shuttle orbiter's return from orbit 
is marked by three mission phases whose flight 
dynamics and constraints are fundamentally dif- 
ferent. (Fig* 1 *) The entry phase occurs at angles 
of attack (o) greater than that required for max- 




figure 1. Flight phases during the return from 
orbit. 


imum lift- drag ratio, L/D, (the "backside of the 
L/D curve") in order to meet heating constraints. 
This phase passes through a region of radio 
frequency "blackout" during which time ground- 
based navigational aids cannot be used; it terminates 
at mesospheric altitude and supersonic velocity. 
Large range and cross-range maneuvers can be 
made during the entry phase. The entry phase is 
matched to the terminal area and landing phase by 
a transition phase, which further decelerates the 
vehicle to subsonic velocity as altitude decreases 
to approximately 40,000 ft. Angle of attack shifts 
from the backside to the frontside of the L/D curve 
during transition, the principal constraints being 
stability and control boundaries, 1 Ranging control 
can be continued during this phase, although it could 
be necessary to limit or suspend' ranging in flight 
regions where attitude control is particularly diffi- 
cult* The shuttle's terminal areaand landing phase 
will be similar to that of standard aircraft opera- 
tions, with the exception that provisions must be 
made for unpowered flight. 

Trajectory and guidance computations for the 
transition phase are complicated by the rapidly 
changing flight environment and by the interdepen- 
dence of state and control* Mach number, angle-of- 
attack, and air-density effects are strong, and their 
ranges of values on any given transition trajectory 
are large* Whether it is used as an in-line part of 
an explicit guidance law or as a source of a priori 
nominal profiles fora simplified guidance scheme, 
afast, reliable means of generating transition flight 
paths and control histories which meet engineering 
and navigational requirements will be useful, 

A numerical approach to defining optimal 
gliding trajectories is presented in this paper* The 
optimization problem ismadeeasier by simplifying 
the model of flight dynamics and by substitution of 
variables* This allows the number of state and 
adjoint variables to be integrated to be reduced. 
Velocity is replaced by kinetic energy, and total 
energy takes the place of time as the independent 
variable. The requirement for descending to a 
particular altitude with a given velocity allows the 
open end-time of the original problem to be replaced 
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by fixed terminal energy, eliminating the additional 
complexity of end-point adjustment. A steepest- 
descent algorithm with near-optimal step-sizing 
rapidly minimizes a cost function consisting of 
integral and terminal penalties. The integral penal- 
ties provide phugoid damping (minimizing dynamic 
pressure- and load factor peaks) and assure that a 
remains within stability and control boundaries* 
The terminal penalties (which are soft constraints) 
force the final state into the desired region. As is 
to be expected with any iterative solution of equa- 
tions, the speed of convergence is greatly affected 
by the choice of initial control profile; however, it 
appears that a small family of starting control 
histories will initiate convergence to an acceptable 
flight path within a few iterations (i.e*, less than 
half a dozen)* 


DEVELOPMENT OF EQUATIONS 

Transformation of Variables 

The equations of motion for the two-dimen- 
sional trajectories considered here make use of the 
flat -earth approximations, on the basis that glide 
range and altitude changes during the transition 
maneuver are small compared to theearth's radius 
and that velocity is decidedly sub-orbital. With the 
further assumption of an air-density profile, P(H), 
which is exponential with altitude, the equations for 
velocitymagnitude (V), flight path angle (Y), altitude 


(H), and range (R) are 

V = -C D ke'^ H V 2 /2 - g sin y (1 

y = C L ke’^ H V/2 - (g/V)cosy (2 

H = V sin y (3 

R = V cos y (4 


where $ is the inverse scale height of the atmo- 
sphere, g is the average gravitational constant in 
the transition altitude interval, and the constant, 
k, combines reference area (S), vehicle mass (m), 
and reference air density (P )as 


m 

Angle of attack (a), is taken as the control 
variable; it enters the equations through the lift- 
and drag coefficients <C L and C-p), which are also 
arbitrary functions of Mach number (M). 

The kinetic energy per unit weight is 


K - V 2 /2g, < 6 

and it is seen that K can replace V in Eq. (1) to 
(4)* Differentiating Eq. (6) with respect to time, 

K * VV/g (7 


and the revised set of equations is 


K = 

-C D ke"' JH (2gK 3 ) 1/ ' 2 

- (2gK) ! / 2 stny 

(3 

y " 

C L ke'^ H (gK/2) 1/2 - 

<g/2K)*/ 2 cos y 

<9 

4 

H = 

{2gK)^ 2 sin y 


(10 

R = 

(2 gK) 1 / 2 cos y 


(11 


Equations (8) to (11) have no explicit depen- 
dence on time; therefore, time can be replaced as 
the independent variable, and the number of equa- 
tions can be reduced by one. Recent papers have 
replaced time by altitude 2 and flight path angle 3 to 
good effect. There is a potential difficulty in 
choosing K, Y, or H to be the independent variable, 
as there is no assurance that these quantities will 
decrease monotonically with time. In particular, 
a phugoid o s c i llati on in ne ar - h or i z ont al flight causes 
the time- rate- of- change of all these variables to 
change sign* As a consequence, it becomes neces- 
sary to transit singular points, and the control 
profile becomes multi-valued. There is a combina- 
tion of K and H which is guaranteed to bemonotonic 
in gliding flight (although not in powered flight 4 ). 
It is, of course, the total energy per unit weight, 
or specific total energy, 

E = K + H, U2 

which is always dissipated by the effect of aerody- 
namic drag. This is easily seen, because 

E = K + H, d3a 

and Eq. (8) and (10) show that 

E = -Cpke'^gK'V^ 2 <l 3b 

Equation (13b) is always less than zero. The 
derivative 

d{ )/dE = [d{ )/dt]/[dE/dt] s { y (14 

is thus well-behaved everywhere on the gliding path, 
and the differential equation for either K or H can 
be eliminated using Eq* (12) and (14), Choosing H 
to be eliminated, the differential equations of motion 
using specific total energy (E) as the independent 
variable become 

K' = 1 + siny/c^ (15 

y' = f-c L + cosy/ p)'2c d k {16 

R' = -cosy/C D M (1? 

where 

M = ke ’ ® H K {18 


2 



and H is determined from Eq* (12)* p can also be 
written as qS/W, i*e*, the dynamic pressure (q) 
divided by wing loading (W/S)* 

Expressing the vehicle dynamics in terms of 
total energy has several advantages in addition to 
eliminating a differential equation* The time to 
execute the transition maneuver is inconsequential; 
given a set of initial conditions on V, V, H, and R, 
the maneuver must be defined to match final condi- 
tions on the same variables, subject to constraints 
which are functions of the state and control, V and 
H are combined in E; hence, the running interval 
[ E , EjO is fixed and well-defined. Control histories 
and gams are then computed as functions of energy 
rather than time-to-go, This introduces a feedback 
relationship between the vehicle's state and control, 
and it eliminates the problem of running out of gains 
and nominal control profile if time-to-go is exceeded 
before the terminal flight condition is reached (for 
information, time-to-go can be calculated using Eq* 
(13b), since dt = dE/E)* A final advantage is that 
the partial derivatives necessary to describe per- 
turbations about a given trajectory are simpler and 
fewer in number (see following sections). 

Equations (6), (12) and (15-17) correspond 
exactly to Eq* (1-4); further simplifications can be 
made for transition trajectories* It is observed 
that flight path angle rates are negligible, i*e*, that 
load factor is close to 1 during the transition* 1 
Assuming that y = 0 leads to V 1 = 0* The algebraic 
solution to Eq* (16) is then 

cost = C l M ( 19 


and the number of equations to be integrated is 
reduced to 2* Alternatively, the flight path angle 
can be assumed small during thetransition, 1 leading 
to a third-order system with cos V = 1 and sin T “ 
7. Both simplifications are examined further, and 
the optimization continues with the latter approxi- 
mation* 

Second- and Third-Order Dynamic Models 

With negligible Y rates, Eq* (19) applies, and 
the dynamical equations (for Y ^ 0) are 

2 2 * 

K' = ITU- C L * #T) / C d M (20 


R' = -C 


L' D 


(21 


The simplicity of the range equation provides 
a well-known result for gliding flight, namely that 
glide range between two energy- states is maximized 
by flying at This is easily seen, since 


R = 



<- c L / c D > dE 


> E f (22 


Theintegral is maximized when the integrand 

is minimized, which occurs at L/ D max * This result 

depends not on small Y but on negligible Y, and it 

allows C t and C~ to be functions of M as well as 
L D . 


fl . In order to study optimal paths, it is necessary 
to find the sensitivity of the state equations to small 
variations in K, R, and a?* Neglecting C T and , 

d M 

the three non- zero partial derivatives are 


bK' „ l 

— -- = + - 


r mC D C L C L " C L ** * C D 


d* (C D \) L (1 - C L V> i'2" 

SR' „ _ / c D C Lg- C L C Dg \ 


d O' 


9K f . ± ^K_ 


SK 

where 


i - C L 2 (1 - a 2 ) "| 

(i - c L W 2 


H K - M ( P + l/K) 


<23 


(24 


(25 


06 


These derivatives are well-posed for small 
perturbations about gliding flight paths, with one 
exception: the denominator of Eq* (23) and Eq. (25) 
can lead to difficulty* The denominator is actually 
sin V; hence the partials of K* become infinite in 
horizontal flight* When C > 1, Eq* (19) is no 
longer valid, and the denominator becomes the 
square root of a negative number* It is likely that 
both difficulties will be encountered in a numerical 
optimization - if not on the optimal path, on a 
previous iteration. For these reasons, the second- 
order model will not be pursued; however, it is 
interesting to note that for planar, gliding trajec- 
tories without integral or terminal- range penalties, 
the equations required for optimization (see next 
section) are scalar* With aterminal-range penalty, 
the range adjoint variable is non-zero but constant, 
leaving a single differential equation to be solved 
for the kinetic energy adjoint variable and a two- 
component control correction. The conventional 
assumption of a "parabolic drag polar" (Ci linear 
in a, C D quadratic in Cl) simplifies Eq, 120-25) 
further, although the assumption is valid only for 
subsonic flight, , 

In the second approximation, the rate-of- 
change of flight path angle is not limited, but y is 
assumed small* Equations (1 5-1 7) can be expressed 
as the vector equation 

xj = _f { x , <y) (27 


where Xj = K, - Y, = R, and 


f ! = 1 + X 2/ C D^ < 28 

f 2 ■ <-C L + 1/4)/2 C d x 1 (29 

f 3 = - 1/Cj^i (30 

Equation (30) showsthat C D Mhedrag decel- 
eration, must be minimized to maximize range. 
Again assuming that and are negligible, 

7 of the 12 partial derivatives with respect to x 
and at are non-zero. They are 
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f l 

= -X2(8 + . 1/x^VCjyj 

-(31 

1 





(32 

f 2 

*1 

= [ C l -(2 + pxj)/#!] /2 C d Xj 2 

(33 

f 3 

X 1 

= ( 6 + I/XjJ/Cjj#! 

(34 

f l 

a 

= " X 2 C D J C D ** 

(35 

f 2 

& 

= [-C, + C D (C L - l/*i)/C D l /2 C dXi 
1-1 a o' 

(36 

f 3 

= c D /C d 2 „ 

(3 7 


The principal virtues of the small-angle ap- 
proximation are computational, in that the number 
of partial derivatives is reduced, and there are no 
trigonometric functions or square roots to be evalu- 
ated* These attributes, plus the coincidence of 
similar terms in Eq, (28-37), make the third-order 
model hardly more complex than the second-order 
model* With range open, the state dimension is 
two, A terminal-range penalty increases the state 
dimension to 3, but the range adjoint variable is 
again constant, leaving two adjoint variables to be 
evaluated by numerical integration. 

Opti mization Using the Third-Order Model 

It is desired to minimize a cost function 
consisting of terminal and integral penalties, sub- 
ject to the dynamic constraint of Eq. (27), The 
augmented cost function, J, for fixed end condition 

J = <x f -x D ) T 9<X f -x D ). 

E f 

+ 5 U(X. #) +X T [f (x.^-x'l } dE, E o >E f (38 
E o 

where Q is a constant, diagonal matrix weighting 
the squared-error between achieved- and desired 
final state, ^ is a penalty function whose integral 
must be minimized, and X is the vector adjoint of 
x* The state vector is a function of E through Eq, 
(27), while ME) is found from 

X'(E) = - f x T (E) X(E) -,>) {3g 


with 


_A(E f ) = 2Q (x f - x D ) (40 


Having obtained a trajectory from Eq* (27) 
with an initial control profile,o?(E)* the angle-of-at- 
tack history is improved on succeeding iterations 
by the perturbation 

flor(E) - - + W 


where € is a near-optimal step-size obtained by a 
one- dimensional search of J(e), Aquadratic approx- 
imation, J(€), is found by evaluating the costs of 
three trajectories with control profiles a (E), a (£) 
+ ^ 0 6 oKE), and 0 ! o (E) + 2 € q 6q:(E), The iteration is 
completed byjchoosing the e which minimizes the 
variation of J with respect to e* Additional test 
trajectories are computed if it is apparent that a 
significant improvement can be made, e,g*,if 

€ (min J) > 2c * 


Although terminal- range corrections are 
most readily made by varying a during the early 
portion of the trajectory, it is found that no combina- 
tion of penalty weights provides this obvious control 
correction. Equation (41) tends to reach its largest 
values in the latter portion of the flight, A simple 
artifice which retains the essential ingredients of 
the steepest-descent algorithm overcomes the 
problem* On the first one or two iterations 
(depending on terminal-range error), the step-size, 
e Q , is multiplied by a ramp function, which equals 
1 at E 0 and 0 at Efj hence* the corrections are 
attenuated astheterminal point is approached. The 
ramp weighting allows large changes in range with 
little change in terminal V and y f which are pri- 
marily determined by the control profile in the 
latter portion of the flight* 

The integral penalty function t£) contains 
terms which enforce control boundaries and which 
introduce trajectory damping* Angle -of- attack 
boundaries imposed by elevator "control power 11 
and embedded regions of instability which are 
functions of Mach number have been considered 
previously* 1 In the present paper, the a 'penalty 
is simply 


1 

,2 

C - Of^) 

. a < 


^1 

{ 0 

/ i v2 

I C 2 ( Of- ar 2 ) 

- a > a 2 

(42 


with of- and Of^ chosen as 3* and 50 * respectively 
(c^ ana C 2 must be negative for the energy integral 
of =£ to be positive). The upper limit roughly 
corresponds to the highest hypersonic trim angle- 
of- attack likely to be obtained for a delta- wing space 
shuttle with negative- maximum elevon deflection. 
The trim limit is sensitive to center-of-gravity 
location and could be sharply reduced by aft move- 
ment of the c,g*; it also decreases with decreasing 
M* The 3*-limit is arbitrarily chosen to maintain 
positive load factor* is academic for the numer- 
ical results which follow. It comes into play during 
numerical iteration, but none of the optimal profiles 
presented here follow a control boundary, 

Phugoid oscillation, which results from a 
lightly damped interchange of kinetic- and potential 
energy, 5 leads to load factor and dynamic pressure 
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overshoots, as well as y and pitch attitude (0) 
oscillation. The oscillation is forced by the continu- 
ally changing flight condition; it can be reduced only 
by a varying a( which may occur somewhat automati- 
cally in planar cases by redefining the control 
variable to be 6 - T-a, as in Kef, 6), Penalties on 
load factor, y, or dynamic pressure can be employed 
to reduce the phugoid oscillation. 

Load factor and y penalties are closely re- 
lated, The load factor is a function of both Ct and 
C D , while (from Eq. (2)), y is a function only of 
C^. Analogously, a penalty on y' 2 using Eq. (29) 
can be formed. The penalty has been used with 
some success, principally because it is a function 
of both the state and the control. 


hypersonic L/D of 2.1, which occurs at & = 13.8 s , 
and a subsonic L/D m g LX of 4.3 (a = 8.4 s ); the wing 
loading is 49 psf. The nominal starting point for 
these trajectories occurs at an altitude of 150,000 
ft and a velocity of 8,000 fps, corresponding to a 
specific total energy of 1.15x10^ ft; the nominal 
terminal point at 40,000 ft altitude and 870 fps 
velocity (M = .9) yields a specific total energy of 
5.18X10 4 ft. At the beginning of the transition 
trajectory, kinetic energy makes the larger contri- 
bution to total energy, whereas potential energy is 
the larger contributor at the end point. Flight is 
nominally horizontal at the start and ends with a 
flight path angle of -18*, which is representative 
of the equilibrium angle during the terminal area 
glide. 


Dynamic pressure (q) is less successful in 
providing trajectory damping, because it is an 
implicit function of the control. The dynamic 
pressure can be expressed as 




(43 


and it is clear that q is an explicit function of the 
state only. Similar difficulties with state constraints 
can be solved by differentiating the constraint with 
respect to the independent variable until the control 
appears. 7 For dynamic pressure, a single differen- 
tiation is sufficient, as 

q' a q <P Xj’ + p' *j) 

= q (( ? Xl + lK f Mjl (44 


and fj is a function of or through C c (Eq, (28)). 
Equation (44) shows thajt the q f penalty is primarily 
a kinetic energy rate (x^) penalty which is weighted 
by air density. The y 1 penalty damps angular 
orientation of the velocity vector, while the q* 
penaltyprovideskineticenergy damping. The latter 
has been found to be more effective, and the trajec- 
tory damping penalty function used here is 


^2 ‘ c 3 q/2, c 3 < (45 

with the partial derivatives 

^2 = 2 c z qq 't * l , x 1 )[(px 1 + 1 ) fj- BxJ 

X 1 

+ [Mf r 1) + Ox + l)f 1 J . (46 

1 

,£ 2 = 2c 3 qq' ( Bx +1) f (47 

APPLICATION TO SPACE SHUTTLE TRANSITION 


The numerical results presented in this sec- 
tion demonstrate a variety of planar transition 
trajectories fora delta-winged configuration of the 
space shuttle orbiter, 8 The vehicle has a maximum 


These results show the detailed effects of 
range control, a comparison of the flat- earth trajec- 
tories with round-earth paths, and variations due 
to changes in initial conditions. A convergence 
example is presented with a discussion of the 
feasibility of "real-time* 1 optimization. It will be 
shown that the "jump" to the "front- side of the L/D 
curve 11 need not be made until the very end of the 
transition, that the dynamic model provides a con- 
servative range estimate, and that the present 
formulation suggests several alternatives for real- 
time trajectory control. 

Trajectories to a Given Range 


This section concentrates on the numerical 
results achieved for five terminal-ranges, which 
are bounded by the L/D^^ case at the upper end 
and by a maximum load factor of 3 11 g T s r \approxi- 
mate)atthe lowerend. Terminal range varies from 
250 nmi to 402 nmi for these cases, all of which 
have the same initial and final specific energy. 
The 350 nmi case closely approximates the range- 
open optimum. It will be seen that the phugoid 
oscillation is not entirely eliminated by the dynamic 
pressure-rate penalty - a point of diminishing 
returns was inevitably reached in the iterative 
process. Because the phugoid is a natural mode 
of response, it is put to good use matching initial 
and final conditions. Furthermore, with the excep- 
tion of the L/ D max case, the maximum dynamic 
pressure peaks never exceed 234 psf, hardly more 
than the 190 psf of the nominal terminal condition 
and not a constraining factor for structure or 
control. 

These trajectories are summarized in Fig. 
2, which presents the variation of velocity with 
altitude. The H-V profile is a direct indication of 
the energy distribution on the transition flight paths. 
While flight at L/D max is generally considered to 
provide a smooth dynamic environment, it evidences 
the largest phugoid oscillations and dynamic over- 
pressures of this section, largely due to the initial 
flight condition. The initial condition is not matched 
to an equilibrium glide, and the vehicle falls into a 
maximum q of 425 psf, A brief pitch-up at the 
beginning of the transition phase absorbs the oscilla- 
tion and provides a flight path with a maximum q 
of. 234 psf, losing 7 nmi of range in the process. 
The 402 nmi range of the L/D max case could not 
be improved upon in the numerical optimization. 
Comparison of the five trajectory profiles with the 
constant-q curves plotted in Fig. 2 illustrates the 
general result that increasing range corresponds 
to increasing q, while the intersections with con- 
stant-E curves show the corresponding result that 
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Figure 2* Altitude- velocity profiles for trajectories 
to several terminal ranges* Initial- and final specific 
energies are fixed for these 5 cases* Phugoid 
oscillations cause dynamic pressure peaks* Range 
increases correspond to higher dynamic pressure* 


(ignoring the phugoid mode) long-range trajectories 
have proportionately higher kinetic energy at a given 
specific energy level. One sees evidence^of the 
"zoom- climb/zoom-dive" and maximum- E seg- 
ments which arise from analytic solutions for sim- 
ilar cases, 2 although the effects of finite y rate-of- 
changeand fixed terminal range cannot be ignored. 
The H-V trajectories converge to similar profiles 
near the end point, suggesting that the optimal 
terminal transition maneuvers are relatively 
independent of range. 

That terminal state and control histories are 
similar is verified by Fig. 3, 4, 6 and 7* TheL/ D max 
case is not constrained to a particular terminal 
flight condition (other than identical specific 
energy); hence the terminal velocity is 50 fps lower 
and the altitude is 1300 ft higher than nominal* 
The maximum terminal variation from nominal for 
the remaining four cases is less than 1/6 of these 
numbers, while the average variation is consider- 
ably less. 

Figure 3, which presents a and $ as functions 
of E, illustrates the large initial separation in angles 
for differing terminal- ranges and the convergence 
on similar profiles as the end point comes near 
Notice in particular that a is virtually always greater 
than the a required for L/Djn a? {Fig. (3a)) and that 
abrupt changes in 0 (Fig. 3©) occur near E = 10 5 
ft, corresponding to an M of about 1,2 and altitudes 
of 75,000 - 80,000 ft* The phugoid oscillation is 
generally apparent in 0 but not n; because 7=0- 
o, it is seen in flight path angle as well. The two 
shortest range cases require high a at the start of 
transition to dissipate specific energy and shorten 
the range. Figure 4(b) shows that maximum load 
factors are associated with the range adjustment 
early in the transition, are close to 1 for a large 
portion of this mission phase, and show an increase 


associated with the terminal pull-up to the desired 
flight path angle (yreaches a minimum of -21,3*). 



Figure 3. Angle of attack (a) and pitch angle (0) 
as functions of the specific energy (E) for several 
terminal ranges* a remains on the "back side of 
the L/D curve" for all cases. 



b) Load Factor 


Figure 4* Dynamic pressure (q) end load factor 
("g's") as functions of E for several terminal 
ranges* The greatest q excursions occur for flight 
at L/ due to non- equilibrium starting 

conditions. 
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Figure 5* The relationship between the logarithm 
of specific energy (E) and time is nearly linear, 
with the slope determined by final range* 


■j 




Figure 6* Lift* drag ratio (L/D) and drag coefficient 
(Cp) as functions of E for several terminal ranges. 
Drag modulation is seen to be the principal means 
of range adjustment. 


The qualitative relationship between specific 
energy and time is simple; the logarithm of E 
decreases nearly linearly with time, and the approx- 
imate slope is a function of the final range. Figure 
5 illustrates the precise relationship, which contains 
some oscillation as well as a small mean curvature. 

Figure 6 shows that the principal means of 
adjusting range is in modulation of the drag coeffi- 
cient (Cp)* The lift-drag ratio changes with Cp 
and is instrumental iii maintaining a smooth flight 
path. To this extent, the shape of the Cp profile 



Figure 7* Angle of attack (or) vs* Mach number (M) 
for several terminal ranges* Stability boundaries 
are typical* 


is determined by L/D modulation; however, the 
average magnitude of Cp is a more significant 
function of desired final range than is L/D* 

A plot of a vs* M (Fig. 7) again shows that 
there is no important requirement for flight on the 
front side of the L/D curve* Flight at lower a 
requires higher q, a condition to be avoided if 
dynamic pressure peaks are a major concern* As 
shown previously, 1 the final reduction in a to the 
trim-glide angle occurs at nearly constant, subsonic 
M, Mach number actually decreases to *85 before 
taking its final value of .9* The n-M plot finds its 
greatest use in defining the stability and control 
environment, as both parameters have strong ef- 
fects on vehicle aerodynamics* Typical boundaries 
along which short-period and Dutch roll modes of 
longitudinal and lateral- directional angular motion 5 
become statically unstable are superimposed on Fig* 
7* Thea-M profiles are shown to avoid longitudinal 
instability, but flight in a region of lateral- direc- 
tional instability is unavoidable for M above 2 to 
3. Stability augmentation is therefore required for 
the lateral- directional mode* 

Errors Due to the Simplified Dynamic Model 

Elimination of the round-earth terms in the 
equations of motion and the small y assumption lead 
to errors in trajectory computation* Ignoring cen- 
trifugal force effectively reduces the lift on the 
vehicle, causing pessimistic estimates of range (too 
short) and dynamic pressure (too high)* Gravita- 
tional acceleration is overestimated at the beginning 
of the trajectory and underestimated at the end; 
however, the altitude span is not large, and this 
simplification has small effect* Range is overesti- 
mated by neglecting the earth's curvature* The V 
assumption leads to an overestimate of gravitational 
effects on K' and Y T , while the range calculation is 
increased byassuminga horizontal velocity vector* 

The over-riding effect is the lack of centrif- 
ugal force* A comparison of the simplified calcula- 
tions with a round-earth model which neglects 
oblateness and rotational effects shows range 
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increases of 3% to 5% for ranges of 250 to 4Q2 nmi 
for the latter, while maximum dynamic pressure 
is reduced from 0% to 16%. Both errors are 
conservative, as the actual trajectory is flown with 
lighter aerodynamic loads, and the terminal-range 
is reached with excess specific energy. Restoring 
centrifugal force to the equations would correct the 
principal error and cause small additional 
complexity, A pragmatic approach would be to bias 
the target range* 

Initial Condition Effects 

Two kinds of initial condition effects are 
considered: those which are accounted for by 

recomputation of the optimal control and those 
which are not* In both cases the control is a function 
of E; hence, there is feedback of altitude and 
velocity, and terminal V, v, and H are close to the 
desired values. The primary difference is whether 
or not range is controlled* This difference is 
reflected in the early control history, and, there- 
fore, in the most significant q peak, 

A variation in initial velocity with range 
controlled to 350 nmi has counter-intuitive effects 
on load factor and dynamic pressure* Variations 
of +500 fps have small effect on maximum load 
factor, but the result is an increase for both positive 
and negative variations* The negative variation 
causes a substantial increase in maximum q, while 
the positive variation shows reduced q* These 
results can be attributed to the phase of the phugoid 
oscillation induced by range control* Lower initial 
velocity leads to lower E 0 , causing the vehicle to 
fly at lower a to stretch the range* As shown by 
Fig, 8, the loss of lift causes H to drop and q to 
rise* The opposite is true when there is excess 
energy. 



VELOCITY, FT.xlO " 3 

Figure 8. The effect of initial velocity on flight 
paths to 350 nmi final range. Low initial velocity 
causes descent to lower altitude, increasing dynamic 
pressure* 


Flight path angle variation (with range con- 
trolled) has an equally interesting effect* Raising 
the initial y provides ballistic lengthening of the 
trajectory. Angle of attack increases to shorten 
the range, and maximum q is reduced* 


Results obtained with initial state variations 
and a single near-optimaia -E profile for 350 nmi 
final range are summarized in Table L The effects 
of the initial variations are computed for the round- 
earth model* The table shows that final V and 7 
are virtually insensitive to the initial variations, 
and Hf sensitivity is small* Scheduling a as a 
function of E (and thus H and V) provides effective 
closed-loop control of the terminal altitude and 
velocity vector. The same cannot be said for final 
range, which shows small changes with initial alti- 
tude variation, moderate change with V, and large 
change with V (which, in this case, represents a 
large change in E 0 ). The increase in maximum 
load factor is greatest with the -T variation, and 
the +y variation provides a 40% increase in maximum 

q. 



v r 

y V 

H f * 

R f , 

* 

Tnax' 

Load 

E , 

o . 

Case 

fps 


n 

nmi 

pgf 

factor 

fUlO 

Nominal (Flat-Earth) 

06 & 

-18.2 

40052 

350. 1 

171 

1.47 

1. 15 

Nominal U*outuI- Earth) 

054 

-16.5 

40458 

366. 

152 

1.46 

1. 16 

+500 fps 

354 

-t6. 5 

40440 

405.6 

120 

1.44 

1.29 

-500 fps 

&54 

-16*5 

40465 

327.9 

157 

1.48 

1.04 

+3° 

054 

-16.5 

40464 

363* 4 

211 

1.59 

1. 16 

-3* 

555 

-16. 7 

40569 

348* 1 

20] 

1*66 

1. 16 

+5000 ft 

053 

-16. 5 

404 66 

360.7 

172 

1.36 

1*16 

-5000 ft 

054 

-16.5 

404 68 

363. 1 

182 

1.64 

1* 15 


* Excluding end points 


Table I* Effects of initial condition variations on 
a 350 nmi transition trajectory* 


Real-Time Applications 

The alternatives available for "real-time 11 
optimal trajectory control fall into three catego- 
ries: on-line solution of the two-point boundary- 

value problem, linearization about a nominal flight 
path to obtain optimal feedback gains* and dynamic 
programming* The present investigation was moti- 
vated by the first alternative but could be useful 
for the remaining two alternatives as well* In any 
case, a practical guidance scheme must include a 
cross- ranging capability, whichisnot studied here. 

Solving the two-point boundary value problem, 
e,g*, the trajectory and optimization equations of- 
fered in this paper, is a formidable task even with 
a ground-based general purpose computer; how- 
ever, the possibility shouldnot be dismissed without 
study, for the flexibility that this approach offers 
is extreme* It is reasonable to assume that only 
one or two full solutions would be required during 
the flight as ameans of improving the nominal flight 
path and control for a perturbation guidance law. 
Nevertheless, integrations along the flight path 
must be fast* and the process must converge reli- 
ably within a few iterations. 

The program which computed this paper's 
numerical results was not "streamlined" for use 
in a flight computer, but it gives some indication 
of what might be expected in a real-time environ- 
ment* Using a variable- energy integration step and 
single-precision arithmetic, about 100 Runge-Kutta 
steps were used for each state and adjoint integra- 
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tion* Near- optimal control-' step adjustment 
required the evaluation of three state histories and 
one adjoint history for each iteration. Fast conver- 
gence is dependent on taking anear-optimal control 
step on each iteration, and the size of the best step 
is not easily predicted from the previous iteration. 
There is no question that starting control profiles, 
control steps, penalty weights, and adjustment rules 
must be chosen carefully for fast, reliable conver- 
gence; however, convergence of the sort illustrated 
in Table II has been obtained for a number of cases. 
The starting profile had V, 7 , and R errors of 35 
fps, 2*7 b , and 16 nmi, with q ma = 181 psf; after 
two iterations, the errors were 2 fps, ,3*, and .04 

nmi, with q ^ =157 psf- 

^max ^ 


Iteration 

Velocity 

cost 

r 

cost 

Range 

cost 

cost 

Total 

cost 


Step 

size 

0 

a* a 

*4 

085* 

a. 5 

1002. 

4643. 

— 

1 

a. 6 

*4 

*01 

3*64 

12.6 

230. 

42* 

2 

,02 

*004 

. 0006 

3*66 

3.60 

. 6 

4* 7 

3 

.004 

. 00 V 

.000* 

3.* 7 

3*08 

*oe 

6*4 


Table 11. Convergence example for 350 nmi final 
range* The inner product, <ba, 6 a>, is a measure 
of the proximity to an extremal path. 


The energy formulation is amenable to neigh- 
boring-extremal linearization, which allows op- 
timal feedback gains to be computed for control 
about a nominal path. 7 The gains for K(E),Y(E), 
and R(E) are obtained by evaluating a 3x3 matrix 
Riccati equation. Although the dimension of this 
equation is not large, there are a large number of 
coefficients to be evaluated, many of which are more 
complex than the coefficients of the original non- 
linear equations; hence, the ''ideal 11 combination of 
computing an optimal trajectory onboard followed 
by calculation of optimal feedback gains appears 
impractical. Pre-computed gains scheduled against 
specific energy would be a better solution, espe- 
cially sincea should remain above thea for L/D max 
during the largest part of the transition. 

A family of optimal transition trajectories 
constitutes an autonomous field of extremals which 
can be used for nonlinear feedback control* The 
theory of dynamic programming 10 proves that there 
is aunique optimal control variable associated with 
each pointin the extremal field, in this case defined 
by K (or V), 7, R and E; hence, orcanbe pre-computed 
as an optimal function of these variables and stored 
within the flight computer, A table-lookup algorithm 
then provides the nonlinear feedback control law. 
The present results suggest that a two - parameter 
table, in which the guidance command, otQ t is a 
function only of E and range -to- go (R^q - Rf - R)* 
may be sufficient* 

Figure 9 illustrates a planar guidance concept 
in which a^(E, Rqq) is supplemented by dynamic 
pressure rate damping to make up for the missing 



Figure 9, Block diagram of dynamic programming 
approach to planar guidance (inner loops of flight 
control system not shown). The guidance command, 
ckq, is a function of specific energy (E) and range-to- 

go * R GO*' 


state variables in the nonlinear function. The 
dynamic pressure rate can be derived from Eq, (44) 
and is, therefore, a function of K and 7, The 
7-dependence can be eliminated by using Eq. (7) 
and (13b) instead of Eq. (28) to find fj, introducing 
a requirement to measure V* 

Figure 10 illustrates (in perspective) the 
guidance function, <*q(E , Rqq ), which is derived from 
the optimal trajectories presented in Fig* 2 to 7. 
Choosing a to lie on the guidance surface results 



Figure 10* Guidance function, Oq (E,R GO ), for the 
planar case* 


in anear-optimal trajectory to the end-point, one 
which yields the desired V, 7, H, and R at the 
beginning of terminal area maneuvering, while mini- 
mizing dynamic pressure rate- of- change. The 
,r shadow M of the surface converges as the final range 
is approached. Ranging control is effectively termi- 
nated in this example several nmi from the end point; 
the surface becomes a line independent of Rgo in 
order to provide a good match to the desired Vf, 
Y f , and H^. 

Table T .II presents some results of applying 
the guidance function of Fig. 10 with several initial 
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guidance function are expected to yield advantages 
for this problem as well. 


Initial 

v r 

v r 

B P 

R f , 

q 

max 

Load 

Conditions 



n 

Ami 

pgf 

factor 

Nominal 

846 

- 15, 9 

40637 

350.0 

135 

1.32 

+500 fps 

847 

-78.0 

40645 

350.0 

775 

1.43 

-500 tps 

045 

-18. 3 

4087: 

349.9 

180 

1.62 

+3° 

B46 

- 16.0 

40860 

350. 0 

137 

1,31 

-3° 

840 

-15, 9 

40653 

350.0 

168 

1.46 

+5000 ft 

840 

-18, 0 

40658 

350,0 

138 

1.32 

-5000 ft 

847 

-16.0 

40653 

350.0 

134 

1.31 


^Excluding end points 


Table III. Trajectories to 350 nmi final range 
using dynamic programming guidance function* 


conditions. No trajectory damping is usedonthese 
round- earth paths to a 350 -nmi final range. End- 
point convergence is seen to be very good, and 
dynamic pressure- and load factor peaks (which 
could be reduced with trajectory damping) are not 
excessive. 

CONCLUSIONS 

The planar gliding path of a space shuttle 
orbiter (or other aircraft) has been shown to be 
readily described as a function of specific energy 
rather than time. Transforming the independent 
variable reduces the order of the trajectory 
problem, and introducing kinetic energy as a state 
variable yields further simplification. The change 
of variables provides a fixed end-point for the 
transition trajectory without restricting the final 
time, a decided asset for flight path optimization. 
As presented, the equations also are applicable to 
terminal area maneuvering and landing approach 
(although only the space shuttle’s transition phase 
has been investigated); the equations could be ex- 
tended to hypersonic entry with little difficulty. 

Trajectory optimization by the steepest- 
descent method has provided typical results for a 
delta- wing orbiter configuration. The problem of 
phugoid oscillation has been reduced by defining a 
trajectory damping integral penalty function which 
is quadratic in the rate-of-change of dynamic pres- 
sure, The trajectories are shown to occur almost 
entirely on the "back side of the L/D curve," 
avoiding the high dynamic pressures associated with 
low C^. 

The choice of specific energy as independent 
variable is particularly fortunate for trajectory 
control, in that optimal control profiles establish 
a nonlinear feedback relationship from altitude and 
velocity to angle of attack. A family of optimal 
trajectories forms a field of extremals suitable for 
control via dynamic programming. The number of 
profiles which are necessary for real-time solution 
is small, because a is a strong function of only two 
variables. Extending the real-time control to 3 -di- 
mensional trajectories introduces a new control 
variable (roll angle) and anew guidance parameter 
(crossrange-to-go, orits equivalent). The specific 
energy formulation and the dynamic programming 
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